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NONLINEAR FLUTTER OF A CIRCULAR CYLINDRICAL SHELL 

IN SUPERSONIC FLOW* 

By David A. Evensen 

NASA Langley Research Center, Hampton, Virginia 
and Mervyn D. Olson 

National Aeronautical Establishment, Ottawa, Canada 

SUMMARY 

A nonlinear analysis is presented for calculating the limiting amplitudes of cylin- 
drical shell flutter by using a four-mode approximation for the shell deflection. The 
aerodynamic pressure is approximated by linear piston theory, and the nonlinearity enters 
the problem through the nonlinear shallow-shell equations for the cylinder. The governing 
equations are reduced to four modal equations by applying Galerkin’s method, and limit 
cycle solutions are obtained by the method of harmonic balance. Stability of the limit 
cycles is investigated numerically by integrating the modal equations on a digital 
computer. 

Two types of limit cycle flutter are obtained: (a) two-mode standing-wave flutter 

and (b) four -mode circumferentially traveling-wave flutter. Under most conditions, the 
two-mode standing-wave flutter becomes unstable and transforms into four-mode 
traveling-wave flutter. The analysis indicates that flutter can occur at aerodynamic 
pressures below the linear flutter boundary. This fact may explain why recent results 
indicate a difference between experiments and linear theory for the flutter of cylindrical 
shells. 


INTRODUCTION 

The self-excited oscillation of thin plates exposed on one side to a parallel super- 
sonic airstream is called panel flutter. For some aerospace applications, the prevention 
of this flutter instability becomes the primary design criterion for the structure. 

Although panel flutter is usually associated with flat or curved plates, a thin-walled 
cylindrical shell can also exhibit this type of instability. Extensive reviews of the panel 

This report, with author credits in reverse order, is also being issued by the 
National Research Council of Canada as Aeronautical Report LR-486, December 1967. 


flutter problem and analyses dealing with the flutter of cylindrical shells are given in 
references 1 to 3. Experimental observations of cylindrical shell flutter are reported 
in references 4 and 5. 

A striking feature of the results of reference 5 is the fact that "almost all the flut- 
ter modes observed in these experiments were of the circumferentially traveling-wave 
type." Such traveling waves are not predicted by linear theory, and it was suggested in 
reference 5 that nonlinearities in the shell were responsible for the phenomenon. This 
suggestion motivated the present study, the purpose of which is to demonstrate that the 
circumferentially traveling-wave flutter can be predicted from a nonlinear analysis. 

Some introductory work on the nonlinear flutter of cylindrical shells appears in 
reference 6, which employed a two-mode approximation to the shell deflection. The pres- 
ent study utilizes four modes, the minimum number required to yield traveling-wave flut- 
ter. In addition, the present analysis retains higher order nonlinear terms which were 
neglected in reference 6. Thus, the present study extends the results of previous work 
to include the possibility of traveling-wave flutter and to examine the influence of higher 
order nonlinearities. 

The aerodynamic pressure acting on the shell is approximated by linear piston 
theory, which is commonly used for the high Mach numbers of interest in this study. At 
the present time, no completely acceptable aerodynamic theory is available for cylindri- 
cal shell flutter. However, in reference 6 it is shown that flutter predictions based on 
piston theory correspond fairly close to experimental findings, at least for moderate 
amounts of internal pressure in the shell. Furthermore, piston theory is by far the 
simplest theory to apply to flutter calculations. Note that there is no real contradiction 
in using linear aerodynamic theory for nonlinear flutter calculations, because the limiting 
amplitudes of flutter will still be within the linear range from the aerodynamic point of 
view. 

The cylindrical shell is represented by the nonlinear shallow-shell equations, and a 
Galerkin procedure is used to reduce the problem to four coupled ordinary nonlinear dif- 
ferential equations for the modal amplitudes. Approximate limit cycle solutions to these 
modal equations are obtained by the method of harmonic balance. These approximate 
limit cycles were verified and their stability was investigated by numerically integrating 
the modal equations on a digital computer. 

The limit cycle solutions for flutter fall into two categories. The first type of flut- 
ter involves only two modes, whereas the second type involves four modes responding 
together in a very special way. The two-mode flutter is similar to the standing-wave 
flutter investigated in reference 6 (and, in fact, is very similar to the flutter of flat 
panels), whereas the four-mode flutter involves circumferentially traveling waves. This 
latter phenomenon is characteristic of axisymmetric structures. A noteworthy feature of 
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the four-mode limit cycle results in the conclusion that flutter can occur below the bound 
ary predicted by linear theory. This observation might explain why cylindrical shell flut 
ter has been found experimentally to occur below the theoretical linear flutter boundary. 

SYMBOLS 


free-stream speed of sound 


A l(t),A 2 (t)\ 

B 1 (t),B 2 (t)j 

a l(t),a 2 (t)l 

bl(t),b 2 (t)j 


generalized coordinates of assumed modes 


nondimensional generalized coordinates; aj = Aj/h, 


A,B,C,D limit cycle amplitudes (see eq. (8)) 

D s shell bending stiffness, Eh 3 / 12(1 - ^2) 


Mo 


Li 


Young's modulus 


N x = 


1 8%F 
R 2 80 2 


stress function, defined by <N X $ = 

N 0 = 


1 9 2 F 

R 9x BQ 

9 2 F 


V. 


9x 2 


aerodynamic pressure parameters (eq. (6)) 


thickness of shell wall 


integers, ranging from 1 to 4 


Kj(xi) functional notation to denote equations (11) or (15) 


ki . 


coefficients defined in appendix A 



L 


length of cylindrical shell 


Moo free-stream Mach number of airstream 

m,n axial and circumferential wave numbers 

N x ,N$ stress resultants due to applied loads 
p pressure 

Poo free-stream static pressure 

p m shell internal pressure (net) 

R radius of cylindrical shell 

t time 

u,v,w median surface displacements of shell (fig. 1) 

x,8,z coordinates in the axial, circumferential, and radial directions, respectively 

(see fig. 1) 

x-l vector representation of the unknowns B, co, f, and cp (see eqs. (C2)) 

x io initial approximation to solution of equations (C2) 

y gas constant, 1.4 

6 Xi corrections to vector xj 0 

A,A',£ damping parameters (see eq. (12)) 

e nonlinearity parameter (eq. (7)) 

K 1 3 L / 8 M oq a or:, 
v Poisson's ratio 


4 



shell material density 


o aspect ratio of shell displacement, 7 rR/nL 

steady-state phase angles 
w,w 0 flutter or vibration frequencies 

C0 ln’ Ct, 2n linear undamped modal frequencies (appendix B) 

A dot over a quantity indicates differentiation with respect to time. 

THEORY 

In the following sections, the governing equations for a cylindrical shell are pre- 
sented and briefly discussed. A four-mode solution is assumed, and Galerkin's method 
is applied to reduce the problem to one involving ordinary differential equations. Approx- 
imate limit cycle solutions to these equations are obtained by the method of harmonic 
balance. Stability of the limit cycle flutter is also discussed. 

Governing Equations 

The nonlinear shallow-shell equations form the starting point for this investigation. 
These equations make use of Donnell’s approximations for cylindrical shells and are 
derived in reference 7. In terms of the radial deflection w and in a common notation, 
these equations are 

9 2 F 9 2 w ,9 2 F8 2 w\ 

at 2 8x2 r2 902 r ax 2 R 2 \8<9 2 8x 2 8x 80 8x80 9x 2 bq 2 J 

( 1 ) 

and 

iv4 F = .i^ + l 

Eh R 9x 2 r2 

where w is the radial deflection of the shell and F is the usual stress function. The 
shell geometry and coordinate system are shown in figure 1. 

The aerodynamic pressure p is approximated by first-order piston theory, which 

yields 


/ 8 2 w \ 2 

8 2 w 8 2 w 

\8x 80/ 

8x 2 80 2 


( 2 ) 
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( 3 ) 



/ 1 9w 3w\ 

P " yP °°^"8F + 00 9 x] 

In writing equations (1) to (3), the reader will 
note that (1) in-plane inertia terms have been 
neglected, and (2) the equations are written for 
w defined positive inward (fig. 1) and p 
defined positive outward. 

Four -Mode Approximation — Galerkin's Method 

Equations (1) to (3) are solved approxi- 
mately by using a four-mode deflection of the 
form: 


w(x,0,t) = |Aq(t)sin + A 2 (t)sin -^r^jcos n 0 + jBi(t)sin + B 2 (t)sin -^^Jsin nd 

2 2 
+ ^jAiWsin ^ + A 2 (t)sin + ^jBj^sin 22 + B 2 (t)sin 


(4) 


These mode shapes are similar to those used previously by the author in the nonlinear 
vibrations of cylindrical shells (ref. 8) and verified experimentally for the nonlinear 
vibrations of rings (ref. 9). The terms which are multiplied by n2/4R are included in 
equation (4) so that the solution will satisfy the periodic continuity condition on the cir- 
cumferential displacement v. 

Substitution of equation (4) into the compatibility equation (2) allows the latter to be 
solved for a particular solution F, which is given in appendix A. The complementary 
solution to equation (2) is taken to be zero. By substituting w and F into the appro- 
priate nonlinear strain-displacement relations of shallow-shell theory, it can be shown 
that the following conditions are satisfied: 

(a) The displacements u, v, and w and their derivatives satisfy periodic con- 
tinuity conditions of the form: 


v(x,0,t) = V(X,0 + 277, t) 

(b) The deflection w goes to zero at the ends of the shell, that is, at x = 0 and 

x = L. 

(c) The boundary conditions for a shell having freely supported ends are satisfied 
to a first approximation. For example, the axial bending moment is composed of terms 
which involve the amplitudes Aj and Bj linearly and nonlinear terms involving 
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products of Ai, Bj, and so forth. At the ends of the shell, the linear terms go to zero, 
but the nonlinear terms do not vanish. A similar situation is true for N x and v, where 
the nonlinear terms do not go to zero at the ends of the shell. Since the nonlinear terms 
are relatively small, it is felt that these boundary conditions approximate those of a 
freely supported shell. 

Finally, the expressions for w, F, and p are substituted into equation (1), and a 
Galerkin procedure is used to obtain four coupled nonlinear ordinary differential equa- 
tions for the modal amplitudes Ai to B 2 . The expressions 9w/9Ai, 9w/9A2, 

9w/9Bj, and 9w/9B2 were used as the weighting functions in the Galerkin procedure. 

In a partly nondime nsionalized form, the resulting coupled equations are 


ai + Aa-^ + (^i n 2 a 1 - fa 2 + -j- a^aiaj + af + bibi + iq 2 ) + -^(a 2 a 2 + a 2 2 + b 2 b 2 + b 2 2 ) 

+ ja 2 + 2a 1 a 2 + a 2 a x + bib 2 + 2b]b 2 + b 2 t>i) + ^^a^aia! + bib^ + 2aj^a 2 a 2 + b 2 b 2 j 

+ 2 a 2 ^a 2 a 1 + a{a 2 + b^ + b]b 2 ^ - ~ 2a 1 ^a 1 a 2 + bib 2 j - a 2 ^! 2 + bj 2 ) + ~ a 2 (a 2 2 + b 2 2 )J 

- ekjaj^aj 2 + b x 2 j + e^ja^aj 2 - bi 2 ) + 4a 1 ^a 2 2 - b 2 2 ) + 2b;^aibi + 4a 2 b 2 jj 

2 r 

- ek 3 a 2 (a 1 a 2 + bib 2 ) + el^a *^ 2 + b 2 2 ) + e 2 i 1 a 1 (a 1 2 + b^) + e 2 i 2 |a 2 (a 1 2 + bi 2 ) 
x (a x a 2 + bib 2 ) + a 1 (a 1 a 2 + bib 2 ) + e 2 ^ 3 a 2 ( a 2 2 + b 2 2 )( a i a 2 + b l b 2 ) 

- 2e ^a^a-j 2 + b^ 2 j^a 2 2 + b 2 2 ^ + e 2 Zga-j^a 2 2 + b 2 2 ^ = 0 (5a) 

a 2 + Aa 2 + w 2n 2 a 2 + fa x + ^- a 2 ja 2 a 2 + a 2 2 + b 2 b 2 + b 2 2 j+ -j 2 jajaj + a! 2 + bibi + bj 2 ) 

+ ^ 1 ^a 1 a 2 + 2a 1 a 2 + a 2 a x + bib 2 + 2^^ + b 2 blj + ^j3a 2 ^a 2 a 2 + b 2 b 2 ) + 2a 2 (a 1 a 1 + bib^ 

+ 2a-j^a;ja 2 + a 2 a^^ + b]b 2 + b 2 b]JJ + ^ a-j^a-^ 2 + bi 2 j + a 2 ^aja 2 + b^b 2 ^ - ^ a i^a 2 2 + b 2 2 j 

+ 4ek 2 j4a 2 ^a 2 2 - b 2 2 ^ + a^aj 2 - bj 2 ) + 2b2(a x bi + 4a 2 b 2 j| - ek 3 a 1 (a 1 a 2 + bjb 2 ) 

(Equation continued on next page) 
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+ ek 4 a 2 ^a 1 2 + t>i 2 j - ek 5 a 2 (a 2 2 + b 2 2 ) + e 2 Z 2 + b l 2 )( a i a 2 + b jb 2 ) + e 2 Z 3 a^a 2 2 + b 2 2 ) 

x ^aja 2 + bjb 2 ^ + a^a^ + bjb^J - e^a^aj 2 + bi 2 ) 2 + 2e 2 Z 5 a 2 (a 1 2 + b! 2 )(a 2 2 + b 2 2 ) 

2 

+ e 2 J 6 a 2 (a 2 2 + b 2 2 ) = 0 (5b) 

bi + Abi + - fb 2 + bi^ajaj + a! 2 + bibi + t>i 2 ) + ^(a 2 a 2 + a 2 2 + b 2 b 2 + b 2 2 ) 

+ -^(a 1 a 2 + 2a.;ia 2 + a 2 a 4 + bjb 2 + 21)^ + b 2 bij + ^^bi^ajaj + bib^ + 2b^a 2 a 2 + b 2 b 2 j 

+ 2b 2 ^aja 2 + a 2 ai + bib 2 + b 2 i>iJJ - ypb^a^ + bjb^ - b 2 ^a x 2 + bj 2 j + 1 b 2 ^a 2 2 + b 2 2 J 

- ekjbj^aj 2 + b^ 2 j + ek 2 Jbi^bi 2 - a] 2 J + 4bj^b 2 2 - a 2 2 j + 2ai^ajbj + 4a 2 b 2 j" 

- ek 3 b 2 ^aia 2 + bxb 2 j + ekjbi^ 2 + b 2 2 ^ + e^jb^aj 2 + bj 2 ^ + e^Jb^aj 2 + bi 2 j 

x + bibrjj + bj^ajag + bib 2 ) 2 J + e 2 Z 3 b 2 (a 2 2 + b 2 2 )(a 1 a 2 + bib 2 ) - 2e 2 Z 4 b 1 (a 1 2 + b x 2 ) 

x (a 2 2 + b 2 2 ) + e 2 Z 5 bi(a 2 2 + b 2 2 ) 2 = 0 (5c) 


\ 6 b 

b 2 + Ab 2 + co 2n 2b 2 + fbj + b 2 ^a 2 a 2 + a 2 2 + b 2 b 2 + b 2 2 ) + + a x 2 + b^ + bj 

+ -^■(* 1*2 + 2a 1 a 2 + a 2 a 4 + bjt^ + 2b 1 b 2 + b^j) + ~|3b 2 ^a 2 a 2 + b 2 b 2 ) + 2b 2 ^a 1 a 1 + bjbjj 

+ 2b 1 ^a 1 a 2 + a 2 a 4 + b^ + b 2 bi)J + y^l^i 2 + b x 2 ) - | b^a 2 2 + b 2 2 ) + 1 b 2 ^a 1 a 2 + bib 2 j 


+ 4ek 2 4b 2 (b 2 2 - a 2 2 ) + b 2 (b 4 2 - a 4 2 ) + 2a 2 ^a 1 b 1 + 4a 2 b 2 J - eksbi^a^ + bib 2 j 


ek 4 b 2 ^a 1 2 + bj 2 ) - ek 5 b 2 (a 2 2 + b 2 2 ) + e^b^aj 2 + bj^^ag + bxb 2 j + e 2 Z 3 

(Equation continued on next page) 
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jbl(a 2 2 + b 2 2 )( a l a 2 + b l b 2) + b 2^ a l a 2 + b l b 2^j - e 2 Z 4 b2^a 1 2 + b] 2 ^ 2 + 2e 2 i 5 b2^a 1 2 + bj 2 ) 

2 

x (a 2 2 + b 2 2 j + e 2 ^6 b 2( a 2 2 + b 2 2 ) = 0 


(5d) 


where a 4 , a 2 , bj, and b 2 are the nondimensional modal amplitudes. The linear 
modal frequencies are designated by and u>2n ( see appendix B), and the aerody- 

namic influence comes through the two parameters 


The nonlinearity parameter e 




A = 


y Poo 

Pghaoo 


f = 


8 y M ooPoc 

3p g Lh 


is given by 


( 6 ) 


e 



(7) 


Note that the equations become linear when e = 0. The coefficients kq to kg and l~i 
to Iq in equations (5) depend upon the shell material properties, mode shapes, and the 
axial membrane stress resultant N x ; they are defined in appendix B. 


Limit Cycle Oscillations — Method of Harmonic Balance 

According to linear flutter theory, disturbances of the modal amplitudes away from 
zero are damped with time as long as the aerodynamic pressure parameter f is below 
the linear flutter boundary, given by f = f 0 . Conversely, when f is greater than f Q , 
the linear theory predicts that disturbances of the modal amplitudes will increase expo- 
nentially with time. However, when nonlinear terms are included in the analysis, the 
situation is altered considerably. The nonlinear terms restrict the growth in amplitude, 
and eventually a steady-state vibration with finite amplitude is usually obtained. Such a 
vibration is termed a "limit cycle oscillation." (See ref. 10 for a discussion of limit 
cycles in nonlinear self-excited systems.) 

Equations (5) exhibit limit cycle oscillations of the type just described. When the 
nonlinearities are relatively small (for example, ekj « w i n ^), the limit cycle vibrations 
are nearly sinusoidal in time. In this case, the method of "harmonic balance” (ref. 11) 
can be used to obtain approximate solutions for the limiting amplitudes and phases. 

To obtain approximate limit cycle solutions to equations (5), the amplitudes a 4 (t) 
to b 2 (t) are assumed to be in the general form: 
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( 8 ) 


aj(t) = A cos wt 
ao(t) = B cos (cot + <p) 

> 

bi(t) = C cos (cot + xf/) 
b 2 (t) = D cos(cot + x) 


that is, the limit cycle vibrations are assumed to be sinusoidal in time. To apply the 
method of harmonic balance, equations (8) are substituted first into equation (5a), and the 
results are grouped into terms multiplying cos cot, terms multiplying sin cot, and higher 
harmonic terms. The terms multiplying sin cot and the terms multiplying cos cot are 
then equated to zero, and the higher harmonics are ignored, since they usually contribute 
little to the solution. 


This procedure results in two coupled algebraic equations involving the unknown 
amplitudes and phases. Applying a similar procedure to equations (5b) to (5d) yields six 
more equations. The end result is a set of eight nonlinear algebraic equations for the 
eight unknowns A, B, C, D, cf> , \f/, x> anc * Particular solutions to these equa- 
tions are now presented. 


Two-mode standing-w ave solution.- A possible approximate solution to equations (5) 
is given by 


aj(t) = A cos wt 
a 2 (t) = B cos(cot + <p) 
bi(t) = 0 

b2(t) = 0 > 


(9) 


which is the same as equations (8) with C = 0 and D = 0. In this case, the nondimen- 
sional shell deflection is given by 


= [A cos cot sin ~ + B cos(cot + 0)sin ^j^Jcos n 6 + -^|a^[ 1 + cos 2cot]sin 2 

+ B 2 Q + cos 2(cot + 0fjsin 2 + 2AB{cos(2ojt + <p) + cos 0jsin — sin (10) 

~ L L L j 


Solutions represented by equations (9) and (10) are referred to as "two-mode standing- 
wave flutter." 

For this type of flutter, the algebraic equations which result from the method of 
harmonic balance reduce to the following four equations involving A, B, (f>, and oo. 
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Au>A + ^ A 2 + B 2 (2 + cos 20)j - |2(k3 - k4 - 4k2) + (^B^sin 2 0^ 

- fB sin 0^ + A 2 + 1 B 2 j + y AB 2 |^ 2 - 4)A 2 + (I3 + ^B^jsin 20 = 0 


(11a) 


(a> 2 - co ln 2 )A + ^ A^3[4(ki - k 2 ) + co^A 2 + 2(2^3 - k4 - 4k 2 ) + o>2)b 2 (2 + cos 2<p) 

+ AcoB^sin 2<^j + fB cos <fi [l + 1 e (| A 2 + i B 2 )J - ^ A ^jA4 + B 2 ^ (l 2 - 4) A 2 
+ (4 + Z5)B 2 J (3+4 cos 20)j = 0 

I b (¥[ a2 < 2 + cos 20) +-^ B^J + |^2 ^kg - k4 - 4k 2 ^ + w 2 jA 2 sin 2<jb| 

Y A 2 b|(4 + £ 5 )B 2 + (l 2 - Z4)A2|sin 20 = 0 (11c) 


(lib) 


AcoB + 


- f A sin 0 


1 + |(| A 2 + -i B 2 
5(4 7 J 


(< 


“ 2 -“2n 2 )B+fj 


B< 2 


1^3 - k4 - 4k 2 j + o) 2 |a 2 (2 + cos 20) + 3^4 (k5 - 16k 2 j+ cu^Jb 2 


Aa>A 2 sin 20> - f A cos 0 


1+ M1 a2+ ^ b2 


- ^ B wZ 6 b4 + A 2 j^2(Z 3 + Z 5 )B 2 


+ (Z 2 - 4)A 2 |(3 + 4 cos 20)^ = 0 


(lid) 


It may be noted that these equations reduce to those already presented in reference 6 
when the e 2 terms are neglected. Solutions to equations (11) were obtained by a gen- 
eralization of Newton’s Method. (See appendix C.) A typical calculation is discussed 
under "Results and Discussion." 

The results were found to be sensitive to small amounts of structural damping. 

The influence of this damping was estimated by assuming viscous damping and modifying 
the parameter A as follows. A new parameter A’ was introduced in place of A, 
where 


A' = A + 2?co ln 


( 12 ) 
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In this expression, A is the aerodynamic damping and £ represents the structural 
modal damping. The new parameter A* was introduced into equations (11) in place of 
A, and the limit cycle calculations were repeated for various values of damping ratio £. 

Under certain conditions, the two-mode limit cycle was found to be unstable with 
respect to disturbances in the quiescent modes b^ and b2. (The procedure used to 
determine stability of the limit cycles is discussed in a subsequent section of this report.) 
For some cases, it was found that if the modes bi and b2 were disturbed from zero, 
their amplitudes increased with time. For nonzero values of bj and b2, all four modes 
participate in the flutter, and "traveling-wave flutter" usually results. Traveling-wave 
flutter is discussed in the following section. 

Four-mode trave ling-wave solution.- A second approximate solution to equations (5) 
is given by 

a^(t) = A cos wt 
ao(t) = B cos(u>t + <fi) 

) (13) 

bi(t) = ±A sin wt [ 

b2(t) = ±B sin(u>t + <p) 

-J 

which is the same as equations (8) with A = C, B = D, \f/ = + and x = 4> + Jr- The 
nondimensional shell deflection corresponding to equations (13) can be written as 

w(x,g,t) _ ^ 7 tx cos ( n g _ ^t) + B sin ^ cos(n0 + <p + u>t) + ■^•jA 2 sin 2 ^ + B 2 sin 2 
h L L 4 L L Li 


+ 2AB sin 22 s in 


^cos,] 


The first two terms of this expression for w/h represent circumferentially traveling 
waves. Hence, this type of flutter is referred to as "traveling-wave flutter." 

For this type of flutter, the algebraic equations obtained by the method of harmonic 
balance again reduce to four equations involving A, B, <p, and co as unknowns. These 
equations are 

AwA - ^03 - 8k2)AB^sin 2<p - fB sin cp 1 - ^r^ 2 - 7 j B 2 ^ + AB 2 |l2A 2 + l3B 2 jsin 20 = 0 

_ — 
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(15b) 


(u>2 - o> ln 2 jA +| A^ki - k2^A 2 + |k 3 - 2k4 + ^k 3 - 8k2^cos 20JB 2 ! + fB cos <p 
X jl +|(a 2 +~ B 2 J| - e 2 A li A 4 - 24A 2 B 2 + b4 +^ 2 A 2 + ~ £ 3 B 2 jB 2 (l + cos 20) 


= 0 


AcuB + ^(k 3 - 8k2^ A 2 B sin 20 - fA sin 0 1+ |^ 2 - ^ B 2 


- A 2 B^ 2 A 2 + £ 3 B 2 )sin 20 = 0 

(15c) 


( 0 * 


X 


" w 2n 2 ) B + | B ^( k 5 " 16k2^B 2 + |k 3 - 2k4 + (k 3 - 8k2)cos 2^^ - f A cos 0 

Z 6 B 4 + 2Z 5 A 2 B 2 - Z 4 A 4 + Z 2 A 2 + Z 3 B 2 |a 2 (1 + cos 20) 


1 +4/a 2 +^ b 2 


= 0 


(15d) 


It is interesting to note that equations (15) are appreciably less complicated than 
equations (11) for the standing-wave flutter. This condition is mainly due to the fact that 
for the assumed traveling -wave solution given by equations (13), all the nonlinear terms 
containing time derivatives in equations (5) cancel out; further discussion of this point is 
given in a subsequent section of the report. 

Equations (15) were also solved by the technique outlined in appendix C; the influence 
of damping was again estimated by using A' in place of A (eq. (12)). 

In addition to the solutions represented by equations (9) and (13), one other type of 
limit cycle flutter was observed. This final type of flutter is designated as "four-mode 
in-phase flutter" and is discussed in the following section. 

Four-m ode in-phase solution .- Another approximate solution to equations (5) is 

a 4 (t) = A cos cot 

a2(t) = B cos(cot + 0) 

> (16) 

bj(t) = ±A cos cut 

b2(t) = ±B cos(cot + 0) 

J 

which is the same as equations (8) with A = C, B = D, 0 = 0 ,n and X = 0 + 0 ,tt. The 
nondime nsional shell deflection corresponding to this solution can be expressed as 
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w(x,0,t) 




7TX 

A cos cot sin — + B cos 
L 


(wt + (p)sin ^j^cos^n# + -^M. 2 (1 + cos 2a>t) 


x sin2 — + B 2 
L 


1 + cos 2(wt + $)Jsin 2 + 2AB|cos(2wt + <p) + cos <^jsin ^ sin (17) 


As may be seen by comparing equation (17) with equation (10), this deflection shape 
is equivalent to that for two-mode standing-wave flutter. The only difference is that one 
deflection shape is rotated circumferentially with respect to the other by +n/4 radians. 
The values of A and B for the two types of flutter are related by the \J2 factor since 
sin(7r/4) = cos(7r/4) = l/j/2. 

The algebraic equations for four -mode in-phase flutter are also related to those for 
the two-mode flutter. One set can be transformed into the other by multiplying the ampli- 
tudes A and B by the \[2 factor. Thus, it can be said that four-mode in-phase flut- 
ter is simply another form of two-mode standing-wave flutter, and the two cases need not 
be treated separately. 

Although the limit cycle solutions obtained by the method of harmonic balance have 
been discussed in some detail, several unanswered questions still remain. First of all, 
there is no guarantee that the three types of limit cycles just described are the only limit 
cycle solutions to equations (5). Possibly other more general limit cycles exist where 
i p is not 0, ±rt/2 or 7 r, where A C, B ^ D, and so forth. Secondly, the question of 
the stability of the various limit cycles has not been investigated. In order to study these 
questions, equations (5) were integrated numerically on a digital computer. 


Solutions by Numerical Integration — Limit Cycle Stability 

The method of Runge-Kutta-Gill was used to integrate equations (5) numerically on 
a digital computer. The procedure used was to begin with a set of initial velocities and 
displacements for a^, a 2 ? bj, and b 2 and integrate forward in time until either the 

amplitudes decayed to zero or some sort of steady-state limit cycle was reached. 

The stability of the limit cycles was also examined numerically. Once a limit cycle 
solution was obtained from the algebraic equations (11) or (15), it was possible to deter- 
mine a set of initial conditions which corresponded to the limit cycle. The numerical 
integration was then started with these particular initial conditions to verify the existence 
of the limit cycle. To test its stability, small perturbations were introduced numerically 
in the initial conditions, and caused them to deviate slightly from those corresponding to 
the limit cycle. Then equations (5) were integrated forward in time by using these per- 
turbed initial conditions. If the numerical integration resulted in the same limit cycle 
solution which was originally perturbed, that solution was said to be stable. However, if 
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the numerical integration did not return to the original limit cycle, then the solution was 
said to be unstable. This procedure was used to determine the stability of the various 
types of flutter motions discussed previously. 

RESULTS AND DISCUSSION 


The preceding calculations were carried out for the shell described in references 5 
and 6. The shell properties, flow conditions, and loading are as follows: 


E = 16 X 106 lb/in 2 (ll X 10 10 N/m 2 ) 
h = 0.0040 inch (0.0001015 meter) 

L = 15.40 inches (0.381 meter) 

R = 8.00 inches (0.203 meter) 
p s = 0.000833 lb-sec 2 /in 4 (8900 kg/m3) 


v = 0.350 
Moo = 3.00 

aoo = 8400 inches/sec (213 m/sec) 

N x = 0 

Ng = 4 lb/inch (701 n/m) 
p m = 0.50 lb/in 2 (3440 n/m 2 ) 

The calculations were nondimensionalized by using the two-mode linear results given by 
equations (C4). Most of these computations were carried out for a circumferential wave 
number n = 23. This value of n was selected since reference 6 shows that it results 
in the minimum static pressure for flutter according to two-mode linear theory. 

Two -Mode Results 

The results of the calculations using equations (11) are shown in figures 2 and 3 for 
the case of no structural damping. Figure 2 is a plot of the variation of the limit cycle 
amplitudes with the static-pressure parameter f/f Q . Note that for small amplitudes A 
and B, the curves of figure 2 approach the linear flutter boundary, which occurs at 
f/f G = 1.0. 
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Modal amplitudes, A,B 



Figure 3.- Flutter frequency and phase angle dependence 
Nondimensional static pressure, f/f 0 on amplitude. Two-mode flutter. Dashed portions of 

curves do not represent real solutions but merely serve 

Figure 2.- Limit cycle amplitudes. Two-mode flutter. to indicate the connection between the two branches, 

n = 23; ? = 0. n = 23; ? = 0. 

In order to interpret figure 2, imagine a shell flutter experiment and consider what 
might happen at various values of static pressure, that is, f/f Q . First, consider 
f/f 0 < 1, that is, conditions below the linear flutter boundary, and imagine that the modes 
aj(t) and a 2 (t) are disturbed from rest. If the disturbances are infinitesimal, they 
will be damped in time and the shell will not flutter. This result is identical with linear 
flutter theory. 

On the other hand, if the disturbances are of the proper magnitude, limit cycle flut- 
ter might occur along segments be (b'c') or de (d’e') of the curves in figure 2. For 
example, imagine that f/f Q = 0.4 and the modes are disturbed such that aj(t) ~ 1.0 and 
a 2 (t) ~ 2.0. Under these conditions, the shell might experience limit cycle flutter at a 
point on the be (b’c') portion of the curves in figure 2. However, numerical stability 
studies have shown that segments be (b'c') represent unstable limit cycles; thus, if this 
limit cycle flutter were disturbed ever so slightly, the amplitudes of aj and a 2 would 
either (a) damp down to zero or (b) increase up to the segments de (d’e’) in figure 2. Seg- 
ments de (d'e T ) represent stable limit cycle flutter, and further disturbances in aj and 
a 2 would just result in a return to the limit cycle, that is, to flutter with A ~ 3.6 and 
B ~ 10.7. 


16 



Now consider what might happen above the linear flutter boundary, that is, at 
f/f Q > 1.0. Under these conditions, infinitesimal disturbances in the modes aj and a2 
increase with time. Linear theory indicates that the flutter amplitudes increase without 
bound; however, with nonlinearities present, the amplitudes eventually stabilize at points 
along ab (a T b f ) or ef (e T f f ) of the curves in figure 2. The branch that will be reached 
(ab or ef) depends upon the initial disturbances as well as f/f 0 . 

It is of interest to note that the lower branches (abc, a T b T c ? ) of the curves shown in 
figure 2 are practically identical with the results given in figure 9 of reference 6. How- 
ever, since terms of order have been retained in the present analysis but were 
neglected in reference 6, the present results have additional branches (def, d T e T f T ) running 
off to the right with increasing aerodynamic pressure. 

In order to define two-mode limit cycle flutter fully, it is necessary to specify the 
flutter frequency and phase angle in addition to the modal amplitudes A and B. Fig- 
ure 3 gives the flutter frequency co/o) 0 and phase angle 0 for various values of the 
modal amplitude A. The curves of figure 3 have been labeled to correspond with 
branches abc and def of figure 2. Note that the dashed portions of the curves in figure 3 
do not represent real solutions; they merely serve to indicate the connection between the 
two branches. The flutter frequency is seen to decrease continuously with increasing 
flutter amplitude, whereas the phase angle first decreases to 90° and then increases at 
the higher amplitudes. Decreasing frequency with amplitude is characteristic of the 
"softening" type of nonlinearity noted in the vibrations of cylindrical shells. (See ref. 8.) 

Figure 4 shows sample time traces of the steady-state motion obtained from the 
numerical integration of equations (5) at f/f 0 = 2.0. To obtain these results, the modal 
amplitudes bi and b2 were set equal to zero in the numerical integration. For 
this particular case, the steady-state motion contains a small amount of beating. The 
influence of the nonlinearities is exhibited by the distinctly nonsinusoidal character of the 
aj curve in figure 4. Despite this nonsinusoidal behavior, the average modal amplitudes, 
the predominant frequency, and the phase angle between the modes obtained from the 
curves of figure 4 agree well with those predicted by the harmonic balance method. (See 
fig. 2, for example.) 

The influence of structural damping on the two-mode results is shown in figure 5. 
For £ = 0 (zero damping), the theory indicates that the limit cycle curves extend to the 
left until f = 0. With damping present, however, the results are fundamentally different. 
The amplitude curves now do not continue to the left until f = 0, but instead bend up 
steeply at small values of f. One would expect that they would then bend again and 
approach the upper branches (def, d T e f f f ) for increasing f. The second-mode amplitude 
B does indeed behave in this way. However, the first-mode amplitude A actually 
overshoots and then approaches the branch def from above. As a consequence, the curves 
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for A exhibit horizontal tangencies at f/f Q = 1.42 and 1.46 for £ = 0.0005 and 0.001, 
respectively. 

Numerical stability studies indicated that the damped curves of figure 5 were 
unstable to the left of the horizontal tangencies just mentioned. For example, with 
£ = 0.001, numerical integration of equations (5) was carried out for f/f Q = 0.38 and 1.20 
in an attempt to obtain limit cycles at these points. Initial conditions for the numerical 
integration were obtained from the curves of figure 5. In both cases, the amplitudes 
decayed to zero and no limit cycle was obtained. On the other hand, when a similar pro- 
cedure was tried at f/f 0 = 1.60, a steady-state limit cycle resulted. These results show 
that the damped two-mode limit cycle solutions have a stability boundary (in terms of 
f/f D ) between 1.2 and 1.6. It is suspected that the boundary occurs at the point of hori- 
zontal tangency (f/fo = 1.46). The boundary location is difficult to pinpoint exactly by 
using numerical integration, since near-zero damping and long integration times are 
involved. 

These stability studies thus identified segments of the damped curves in figure 5 
(to the right of f/f 0 = 1.42 and 1.46) which represent limit cycles which are stable with 
respect to disturbances in the modes a^ and & 2 . When other disturbances were con- 
sidered, however, it was found that the two-mode limit cycles were unstable with respect 
to the modes bi and b2- This instability with respect to the b-modes usually leads to 
traveling-wave flutter. 


Traveling-Wave Results 

Results for traveling-wave flutter were calculated from equations (15) and are 
shown in figures 6 to 8. Figure 6 shows the variation of modal amplitude with 
aerodynamic-pressure curves. These curves are fundamentally different from the cor- 
responding ones for standing-wave flutter in that they first proceed upward with negative 
slope until a vertical tangency occurs at a finite value of the aerodynamic pressure. The 
curves then bend back to positive slope and proceed indefinitely to the right, the first- 
mode amplitude being larger than the second, that is, A > B. Note that B was larger 
than A for the standing-wave flutter at large amplitudes. 

Figure 7 shows how the flutter frequency and phase angle between the modes vary 
with the traveling-wave flutter amplitude. In this case, the flutter frequency first 
decreases very slightly and then increases with amplitude. This result is in contrast to 
the standing-wave flutter, where the frequency decreased with increasing amplitude. Note 
that the amount of frequency shift is very much smaller here than for the standing-wave 
case. The large softening nonlinearity exhibited by the standing-wave flutter may be 
attributed to the nonlinear terms involving time derivatives in equations (5). As noted 
previously, these terms cancel out for the traveling-wave solution, nonlinearities being 
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Figure 6.- Limit cycle amplitudes. Traveling-wave flutter, 
n = 23 : Z = 0. 


left only in the spring and aerodynamic terms 
The final hardening nonlinearity shown by the 
frequency curve in figure 7 results from the 
spring terms in equations (5). The 
phase angle between modes shown in fig- 
ure 7 is seen to decrease at high amplitudes, 
but it does not approach 90° as did the phase 
for the standing-wave flutter. 

The influence of structural damping on 
the limit cycles for traveling-wave flutter is 
shown in figure 8. For clarity, only the 
first-mode amplitude is plotted, since the 
effect on the second-mode amplitude was 
very similar. It is seen that the damping 
shifts the point of vertical tangency slightly 
to the right, but does not alter the nature of 
the curves. This result is in contrast to the 
standing-wave case where a small amount of 
structural damping caused a large change in 
the shape of the limit cycle curves. 



90 ioo "" ' lib’ 120 130 


Phase angle between modes, <£, deg 

Figure 7.- Flutter frequency and phase angle 
dependence on amplitude. Traveling- wave flutter, 
n = 23; £ = 0. 



Figure 8.- Influence of damping on the traveling-wave results, 
n = 23. 
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The results of a numerical integration involving all four modes are shown in fig- 
ure 9 for f/f Q = 2.0. The initial conditions for this calculation were taken from the 
limit cycle results of figures 6 and 7. In contrast with standing-wave flutter, these time 
traces are remarkably close to pure sinusoidal waves of constant amplitude and fre- 
quency. The numerical integration results thus verify the amplitudes, frequency, and 
phase angle predicted for this case by the method of harmonic balance. 




Time, sec 


Figure 9.- Numerical integration results. Traveling- wave flutter; f/f 0 = 2.0. n = 23 ; £ - 0. 


Stability studies involving perturbations in all four modes showed that only those 
portions of the curves in figure 6 with positive slope represent stable limit cycle oscilla- 
tions. The same results would be expected to hold for the cases with structural damping 
(fig. 8). 

Although most of the calculations were made for the n = 23 mode(s), qualitatively 
similar results were obtained for other values of the circumferential wave number n. 
This fact was demonstrated in reference 6 for two-mode flutter, and it is illustrated in 
figure 10 for traveling- wave flutter. The curves of figure 10 were calculated for values 
of n ranging from 17 to 27 and for zero structural damping. Figure 10 shows that the 
flutter mode depends upon the spatial distribution of the initial disturbances as well as 
upon their magnitude. For example, at f = 1.0 x 10^, it is possible to excite limit cycle 
flutter for various n-values, depending upon the circumferential wave number corre- 
sponding to the initial disturbance. If the only disturbances present are infinitesimal, 
figure 10 shows that the shell will not flutter until f = 1.405 x 10^, which is the linear 
flutter boundary for the n = 23 mode. Note that regardless of the n-value excited, the 
minimum limiting amplitudes of flutter given in figure 10 are about 4 to 5 shell 
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thicknesses. This analytical result can be compared with the experimentally observed 
values of 1 to 2 shell thicknesses reported in reference 5. 



Aerodynamic pressure parameter, f 

Figure 10.- Limit cycle amplitudes for various values of n. Travel ing-wave flutter. £ = 0. 

Once limit cycle flutter is initiated for some value of the circumferential wave 
number, it is not clear whether modes with other n-values will begin to participate in the 
motion or not. The experimental results (refs. 4 and 5) suggest that flutter occurs with 
only single values of n. However, the flutter mode (that is, the value of n) may change 
as the static pressure is varied. For example, figure 8 of reference 5 shows that the 
experimental flutter amplitudes do not increase monotonically with increasing stagnation 
pressure; instead, they exhibit noticeable fluctuations. With respect to figure 10, these 
experimental results suggest the possibility that the shell might flutter in one mode 
(n = 22, for example) and then change to another mode (perhaps, n = 23) as the param- 
eter f increases. Such a change in mode might be associated with a jump downward in 
amplitude from the n = 22 curve to the n = 23 curve in figure 10. Additional mode 
changes might occur as f is further increased; however, a thorough discussion of this 
problem would necessitate an analytical study involving combinations of circumferential 
harmonics n. 


CONCLUDING REMARKS 

A four-mode nonlinear analysis of cylindrical shell flutter using piston theory and 
the nonlinear shallow-shell equations has been presented. The results obtained from this 
analysis lead to the following conclusions: 
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1. The occurrence of circumferentially traveling-wave flutter, which has been 
observed experimentally, is predictable by the nonlinear theory. 

2. Limit cycle flutter in a circumferentially traveling-wave mode can occur below 
the stability boundary predicted by linear theory. This fact may help to explain the dif- 
ferences between experimental results and linear theory. 

3. The two-mode standing-wave limit cycle results were found to be very sensitive 
to small amounts of structural damping, whereas the traveling-wave results were not. 

4. The limit cycle amplitudes predicted by the present analysis are somewhat 
larger than the experimental results. It is expected that the theoretical limiting ampli- 
tudes would be reduced if the longitudinal displacement u was made to vanish at the 
ends of the shell. Such an axial restraint would be more representative of the experi- 
ments than are the boundary conditions satisfied in the present analysis. 

5. One limitation of the present study is the fact that it has been restricted to only 
four modes. It is expected that adding more modes in the axial direction would alter the 
present results quantitatively but not qualitatively. On the other hand, it is harder to 
estimate the influence of adding more circumferential modes to the analysis. When the 
aerodynamic pressure exceeds the linear flutter boundaries of two or more modes with 
different values of circumferential wave number n, it is not clear how the modes might 
combine. The answer to this question is left to a future investigation. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., October 31, 1967, 

124-08-05-08-23. 
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APPENDIX A 

PARTICULAR SOLUTION TO EQUATION (2) 


Define a = tt/L and /3 = n/R. Then a particular solution to equation (2) which 
corresponds with the assumed deflection shape w is 



where Aj., A 2 , Bj, and B 2 are functions of time. 
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APPENDIX B 

DEFINITION OF CONSTANTS USED IN EQUATIONS (5) 

The linear undamped modal frequencies are given by 
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Similarly, the coefficients l j to Iq are given by 
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APPENDIX C 


SOLUTION OF THE ALGEBRAIC EQUATIONS 

The complexity of the algebraic equations (11) and (15) precludes the possibility of 
a closed-form solution, and a numerical method must be employed. Such a method was 
introduced in reference 12 and yielded satisfactory results. The method as used herein 
is outlined in the following paragraphs. 

The four unknowns A, B, cp, and o> in equations (11) and (15) are required as 
functions of the aerodynamic pressure for a given shell geometry, speed of sound, 
and Mach number of the airstream, and for various values of the circumferential wave 
number n. The aerodynamic pressure comes into the algebraic equations through the 
two parameters A and f. The damping parameter A* (see eq. (12)) may be put in 
the form: 


A' = A + 2?w ln = kS + 2?w ln (Cl) 

where k = SL/SMcoa^. Hence, the parameter f may be used to represent the aerody- 
namic pressure. 

It is convenient for the calculations required herein to treat A as known and B, 
ct>, f, and (p as unknowns. Let the four component vector (i = 1 to 4) denote the 
unknowns such that x^ = B, X 2 = c o, X 3 = f, and X 4 = (p ; then each set of equations (11) 
or (15) is of the functional form 


Kj(xi) = 0 


(i,j = 1,2, 3, 4) (C2) 


The solutions to these equations are obtained by a generalization of Newton's 
method to many variables as follows. If Xj 0 is a good initial approximation to the 
sought-for solution Xj such that x ^ 0 = xi - Sxj, then by expanding equations (B2) in a 
Taylor series from Xj Q to Xj, 



(C3) 


Now the 6 xj may be found from equation (C3), since by definition Kj(xj 0 ) / 0. When 
this is done, the new approximation for Xj is Xjj = Xj 0 + 6 x^, and the foregoing pro- 
cedure is repeated over again until the corrections Sxj become negligibly small. 
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The initial approximation used to start this process is that obtained from the 
linear flutter problem. This approximation results from the linearization of equa- 
tions (11) or (15) and may be shown to be 


A = B 


~ „ w 1r 2 + w 9 2 

w 2 = o, o 2 = 2n _ 
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Initially, A is taken to be a small number, and the linear solution given in equation (C3) 
is used to start the iteration process. Thereafter, A is increased by a small amount, 
and the solution for the previous value of A is used to restart the iterations. 
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